This script plots individual antibody trajectories in the Asembo, Kenya cohort from enrollment to follow-up. 199 children were measured at both timepoints.
#-----------------------------
# preamble
#-----------------------------
library(here)## here() starts at /Users/benarnold/enterics-seroepi
here()## [1] "/Users/benarnold/enterics-seroepi"
# load packages
library(tidyverse)## ── Attaching packages ─────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.0.0 ✔ purrr 0.2.4
## ✔ tibble 1.4.2 ✔ dplyr 0.7.4
## ✔ tidyr 0.8.0 ✔ stringr 1.3.1
## ✔ readr 1.1.1 ✔ forcats 0.3.0
## ── Conflicts ────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
# bright color blind palette: https://personal.sron.nl/~pault/
cblack <- "#000004FF"
cblue <- "#3366AA"
cteal <- "#11AA99"
cgreen <- "#66AA55"
cchartr <- "#CCCC55"
cmagent <- "#992288"
cred <- "#EE3333"
corange <- "#EEA722"
cyellow <- "#FFEE33"
cgrey <- "#777777"#-----------------------------
# load the formatted data
# created with
# asembo-enteric-ab-data-format.Rmd -->
# asembo-enteric-ab-distributions.Rmd
#-----------------------------
dl <- readRDS(here("data","asembo_analysis2.rds"))Identify children whose status changed between enrollment and follow-up. Those who changed from negative to positive are seroconverters (seroi below), and those who changed from positive to negative are seroreverters (seror below).
#-----------------------------
# identify incident
# seroconversions and reversions
# based on crossing the
# seropositivity cutoff
# and
# based on a 4-fold change in MFI
# to above the cutoff (seroconversion)
# or starting above cutoff (seroreversion)
#-----------------------------
dlp <- dl %>%
group_by(antigen,antigenf,childid) %>%
mutate(seroposA=ifelse(time=="A",seropos,NA),
seroposA=max(seroposA,na.rm=TRUE),
seroposB=ifelse(time=="B",seropos,NA),
seroposB=max(seroposB,na.rm=TRUE),
seroi=ifelse(seroposB==1 & seroposA==0,1,0),
seror=ifelse(seroposB==0 & seroposA==1,1,0),
mfiA=ifelse(time=="A",logmfi,NA),
mfiA=max(mfiA,na.rm=TRUE),
mfiB=ifelse(time=="B",logmfi,NA),
mfiB=max(mfiB,na.rm=TRUE),
seroi4fold=ifelse(logmfidiff>log10(4) & mfiB>serocut,1,0),
seror4fold=ifelse(logmfidiff< - log10(4) & mfiA>serocut,1,0)
) Longitudinal antibody trajectories of individual children for all antigens measured, colored by the change in antibody levels (increases and decreases)
#-----------------------------
# Plot individual antibody
# trajectories, colored
# by increases and decreases
#-----------------------------
# create an change category factor for plot aesthetics
dlp <- dlp %>%
mutate(mficat=factor(ifelse(logmfidiff>0,"Increase","Decrease"),levels=c("Increase","Decrease")))
# convert time to numeric (for plotting)
dlp <- dlp %>%
mutate(svy=ifelse(time=="A",0,1))
log10labs <- c(
expression(10^0),
expression(10^1),
expression(10^2),
expression(10^3),
expression(10^4)
)
ppq <- ggplot(data=filter(dlp,!is.na(mficat)),aes(x=svy,y=logmfi,group=factor(childid),color=mficat)) +
# geom_point(alpha=0.2) +
geom_line(alpha=0.2) +
# geom_hline(aes(yintercept=mixcut),linetype="dashed") +
scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
scale_y_continuous(breaks=0:4,labels=log10labs) +
scale_color_manual(values=c(corange,cteal,"gray80"),
guide=guide_legend(title="MFI change:",override.aes = list(alpha=1))
) +
facet_wrap(~antigenf,nrow=6,ncol=2) +
labs(x="Measurement",y="Luminex response (MFI-bg)") +
theme_minimal() +
theme(
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
axis.ticks.y=element_line(color="gray40"),
legend.position = "top"
)
ppqCreate the same figure of individual antibody trajectories between measurements, but limit it to antigens with seroincidence measures, and color the trajectories by seroconversion/reversion status.
#-----------------------------
# pull the incidence info
# and just merge it back
# to the full longidudinal data
#-----------------------------
dl2 <- dlp
# create an incidence category factor for plot aesthetics
# for seroconversion based on seropositivity only
dl2$icat <- factor(NA,levels=c("Seroconversion","Seroreversion","No change"))
dl2$icat[dl2$seroi==1] <- "Seroconversion"
dl2$icat[dl2$seror==1] <- "Seroreversion"
dl2$icat[dl2$seroi==0 & dl2$seror==0] <- "No change"
# create an incidence category factor for plot aesthetics
# for seroconversion based on 4-fold changes in MFI
dl2$i4cat <- factor(NA,levels=c("Seroconversion","Seroreversion","No change"))
dl2$i4cat[dl2$seroi4fold==1] <- "Seroconversion"
dl2$i4cat[dl2$seror4fold==1] <- "Seroreversion"
dl2$i4cat[dl2$seroi4fold==0 & dl2$seror==0] <- "No change"
# create an incidence category factor for plot aesthetics
# for seroconversion comparing both methods
dl2$icat_comp <- factor(NA,levels=c(">4-fold increase, across cutoff",">4-fold decrease",">4-fold increase, above cutoff","<4-fold change"))
dl2$icat_comp[dl2$seroi4fold==1 & dl2$seroi==1] <- ">4-fold increase, across cutoff"
dl2$icat_comp[dl2$seror4fold==1] <- ">4-fold decrease"
dl2$icat_comp[dl2$seroi4fold==1 & dl2$seroi==0] <- ">4-fold increase, above cutoff"
dl2$icat_comp[is.na(dl2$icat_comp)] <- "<4-fold change"
pp <- ggplot(data=filter(dl2,!is.na(icat)),aes(x=svy,y=logmfi,group=factor(childid),color=icat_comp)) +
# geom_point(alpha=0.2) +
geom_hline(aes(yintercept=serocut),linetype="dashed",size=0.3) +
geom_line(alpha=0.2) +
scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
scale_y_continuous(breaks=0:4,labels=log10labs) +
scale_color_manual(values=c(corange,cteal,cmagent,"gray80"),
guide=guide_legend(title="MFI change:", override.aes = list(alpha=1), ncol=2,nrow=2),
) +
facet_wrap(~antigenf,nrow=6,ncol=2) +
labs(x="Measurement",y="Luminex response (MFI-bg)") +
theme_minimal() +
theme(
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
axis.ticks.y=element_line(color="gray40"),
legend.position = "top"
)
pp#-----------------------------
# final figure, excluding cholera
# due to cross-reactivity with ETEC
#-----------------------------
relev <- levels(dl2$antigenf)[c(1:6,11,7,9,10)]
dl3 <- dl2 %>%
filter(!is.na(icat) & antigen!="cholera") %>%
ungroup() %>%
mutate(antigenf=factor(antigenf,levels=relev))
pp2 <- ggplot(data=dl3,aes(x=svy,y=logmfi,group=factor(childid),color=icat_comp)) +
# geom_point(alpha=0.2) +
geom_hline(aes(yintercept=serocut),linetype="dashed",size=0.3) +
geom_line(alpha=0.2) +
scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
scale_y_continuous(breaks=0:4,labels=log10labs) +
scale_color_manual(values=c(corange,cteal,cmagent,"gray80"),
guide=guide_legend(title="MFI change:", override.aes = list(alpha=1), ncol=2,nrow=2),
) +
facet_wrap(~antigenf,nrow=6,ncol=2) +
labs(x="Measurement",y="Luminex response (MFI-bg)") +
theme_minimal() +
theme(
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
axis.ticks.y=element_line(color="gray40"),
legend.position = "top"
)
pp2# save PDF and TIFF versions
ggsave(here("figs","Fig5-asembo-ab-trajectories.pdf"),plot=pp2,device=cairo_pdf,width=6,height=10)
ggsave(here("figs","Fig5-asembo-ab-trajectories.TIFF"),plot=pp2,device="tiff",width=6,height=10)sessionInfo()## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 forcats_0.3.0 stringr_1.3.1 dplyr_0.7.4
## [5] purrr_0.2.4 readr_1.1.1 tidyr_0.8.0 tibble_1.4.2
## [9] ggplot2_3.0.0 tidyverse_1.2.1 here_0.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.18 cellranger_1.1.0 pillar_1.2.2 compiler_3.5.0
## [5] plyr_1.8.4 bindr_0.1.1 tools_3.5.0 digest_0.6.16
## [9] lubridate_1.7.4 jsonlite_1.5 evaluate_0.10.1 nlme_3.1-137
## [13] gtable_0.2.0 lattice_0.20-35 pkgconfig_2.0.2 rlang_0.2.2
## [17] psych_1.8.4 cli_1.0.0 rstudioapi_0.7 yaml_2.1.19
## [21] parallel_3.5.0 haven_1.1.1 withr_2.1.2 xml2_1.2.0
## [25] httr_1.3.1 knitr_1.20 hms_0.4.2 rprojroot_1.3-2
## [29] grid_3.5.0 glue_1.2.0 R6_2.2.2 readxl_1.1.0
## [33] foreign_0.8-70 rmarkdown_1.9 modelr_0.1.2 reshape2_1.4.3
## [37] magrittr_1.5 backports_1.1.2 scales_1.0.0 htmltools_0.3.6
## [41] rvest_0.3.2 assertthat_0.2.0 mnormt_1.5-5 colorspace_1.3-2
## [45] stringi_1.2.2 lazyeval_0.2.1 munsell_0.5.0 broom_0.4.4
## [49] crayon_1.3.4